      SUBROUTINE RK78(DER,T,TOUT,N,X,H,RELERR,ABSERR,SP)  
C  VERSION OF 4/1/85
C  PURPOSE
C    RUNGE-KUTTA 7(8) METHOD AS GIVEN BY THE FEHLBERG COEFFICIENTS
C    NASA TR R-287
C  INPUT
C    DER    = NAME OF SUBROUTINE CONTAINING THE DERIVATIVES
C    T      = INITIAL INTEGRATION TIME
C    TOUT   = FINAL INTEGRATION TIME
C    N      = NUMBER OF EQUATIONS TO BE INTEGRATED
C    X      = THE STATE TO BE INTEGRATED, DIMENSION 6 HERE 
C    H      = INITIAL GUESS OF STEP-SIZE
C    RELERR = RELATIVE TOLERANCE
C    ABSERR = ABSOLUTE TOLERANCE
C    SP     = APPROXIMATE TIME THAT USEROP WILL BE CALLED
C  OUTPUT
C    T      = AT TOUT, T=TOUT
C    X      = STATE AT TOUT
C  CALL SUBROUTINES
C    DER, USEROP
C  REFERENCES
C    JPL EM 312/87-153, 20 APRIL 1987
C    NASA TR R-287, OCTOBER 1968
C  ANALYSIS
C    J. H. KWOK - JPL
C  PROGRAMMER
C    J. H. KWOK - JPL
C  PROGRAM MODIFICATIONS
C
C  DECLARED DER AS EXTERNAL TO ELIMINATE COMPILER WARNING.
C  F. S. PATT, JANUARY 13, 1995.
C
C    NONE
C  COMMENTS
C    USER MUST CALL RK78CN BEFORE CALLING RK78.
C    THE DIMENSIONS OF XDUM AND F CAN BE CHANGED FROM 6 TO ACCOMMODATE
C    ANY SYSTEM OF DIFFERENTIAL EQUATIONS.  USEROP IS A USER OPTION
C    ROUTINE WHERE THE HEALTH OF THE INTEGRATION CAN BE MONITORED BY
C    SETTING SP=ZERO.  IT CAN ALSO BE USED TO CHECK CERTAIN CONDITIONS
C    OF THE STATE, AND POSSIBLY RESTART WITH A NEW STATE OR STEP-SIZE
C    BY SETTING ISTART=1 IN USEROP.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/FELCON/CH(13),ALPH(13),BETA(13,12),ORDRCP,NORDER,NTIMES
      DIMENSION XDUM(6),F(6,13),X(*)
      EXTERNAL DER
      DATA ZERO,ONE/0.D0,1.D0/
C
C *** FACT IS A SCALING FACTOR TO REDUCE THE ESTIMATED STEPSIZE TO AVOID
C *** EXCEESIVE NUMBER OF REJECTION.
C *** SMALL IS A SMALL NUMBER COMPARED TO THE SIZE OF T
C
      DATA FACT,SMALL/0.7D0,1.D-8/
C
C *** INITIALIZATION
C
      NSTP=0
      NREJ=0
      SDT=DSIGN(ONE,TOUT-T)   
      DT=DABS(H)*SDT
      SPP=DABS(SP)*SDT
      TPR=T
      ISTART=0
      IFLAG=0
C
C *** CALL USEROP AT INITIAL T
C
C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
C
C *** START INTEGRATION AND SAVE INITIAL INFO
C
  100 CONTINUE
      ISTART=0
      NSTP=NSTP+1
      TDUM=T
      TPR=TPR+SPP
      DO 110 I=1,N   
  110 XDUM(I)=X(I)
C
C *** CHECK FOR T+DT PASSING TOUT
C
      IF (DABS(DT).GT.DABS(TOUT-T)) DT=TOUT-T  
C
C *** CHECK FOR REACHING TOUT
C
      IF (DABS(T-TOUT).LT.SMALL) GO TO 900
C
C *** START FUNCTION EVALUATIONS
C
      CALL DER(T,X,F(1,1))
      DO 300 K=2,NTIMES  
      KK=K-1
      DO 310 I=1,N   
      TEMP=ZERO
      DO 320 J=1,KK  
  320 TEMP=TEMP+BETA(K,J)*F(I,J)
  310 X(I)=XDUM(I)+DT*TEMP  
      T=TDUM+ALPH(K)*DT 
      CALL DER(T,X,F(1,K))  
  300 CONTINUE
C
C *** COMPUTE NEW STATE
C
      DO 400 I=1,N  
      TEMP=ZERO 
      DO 410 L=1,NTIMES 
  410 TEMP=TEMP+CH(L)*F(I,L)
  400 X(I)=XDUM(I)+DT*TEMP
C
C *** START LOCAL TRUNCATION ERROR COMPUTATION
C
      ERR=RELERR
      DO 500 I=1,N  
      TER=DABS((F(I,1)+F(I,11)-F(I,12)-F(I,13))*CH(12)*DT)
      TOL=DABS(X(I))*RELERR+ABSERR
      CONST=TER/TOL
      IF (CONST.GT.ERR) ERR=CONST
  500 CONTINUE  
C
C *** ESTIMATE NEW STEP-SIZE
C
      DT=FACT*DT*(ONE/ERR)**ORDRCP
C
C *** IF IFLAG.EQ.1, ACCEPT THIS STEP UNCONDITIONALLY
C
      IF (IFLAG.EQ.1) THEN
        IFLAG=0
        GO TO 510
      ENDIF
C
C *** REJECT LAST STEP IF ER.LT.ONE OR ISTART.EQ.1
C
  220 CONTINUE
      IF (ERR.GT.ONE.OR.ISTART.EQ.1) THEN
        T=TDUM
        DO 610 I=1,N
  610   X(I)=XDUM(I)
        NREJ=NREJ+1
        NSTP=NSTP-1
        GO TO 100
      ENDIF
C
C *** CHECK FOR TIME TO CALL USEROP
C
  510 CONTINUE
      IF (SDT.GT.ZERO.AND.T.GE.TPR) GO TO 200
      IF (SDT.LT.ZERO.AND.T.LE.TPR) GO TO 200
      GO TO 100
  200 CONTINUE
C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
      IF (ISTART.EQ.1) GO TO 220
      GO TO 100
  900 CONTINUE
C
C *** CALL USEROP THE LAST TIME
C
C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
      RETURN
      END
